Phystech@DataScience¶

Домашнее задание 3¶

Правила, прочитайте внимательно:

  • Выполненную работу нужно отправить телеграм-боту @thetahat_ds25_bot. Для начала работы с ботом каждый раз отправляйте /start. Дождитесь подтверждения от бота, что он принял файл. Если подтверждения нет, то что-то не так. Работы, присланные иным способом, не принимаются.
  • Дедлайн см. в боте. После дедлайна работы не принимаются кроме случаев наличия уважительной причины.
  • Прислать нужно ноутбук в формате ipynb. Если вы строите интерактивные графики, их стоит прислать в формате html.
  • Следите за размером файлов. Бот не может принимать файлы весом более 20 Мб. Если файл получается больше, заранее разделите его на несколько.
  • Выполнять задание необходимо полностью самостоятельно. При обнаружении списывания всем участникам списывания дается штраф -2 балла к итоговой оценке за семестр.
  • Решения, размещенные на каких-либо интернет-ресурсах, не принимаются. Кроме того, публикация решения в открытом доступе может быть приравнена к предоставлении возможности списать.
  • Обратите внимание на правила использования ИИ-инструментов при решении домашнего задания.
  • Код из рассказанных на занятиях ноутбуков можно использовать без ограничений.
  • Для выполнения задания используйте этот ноутбук в качестве основы, ничего не удаляя из него. Можно добавлять необходимое количество ячеек.
  • Комментарии к решению пишите в markdown-ячейках.
  • Выполнение задания (ход решения, выводы и пр.) должно быть осуществлено на русском языке.
  • Решение проверяется системой ИИ-проверки No description has been provided for this image ThetaGrader. Результат проверки валидируется и исправляется человеком, после чего комментарии отправляются студентам.
  • Если код будет не понятен проверяющему, оценка может быть снижена.
  • Никакой код из данного задания при проверке запускаться не будет. Если код студента не выполнен, недописан и т.д., то он не оценивается.

Важно!!! Правила заполнения ноутбука:

  • Запрещается удалять имеющиеся в ноутбуке ячейки, менять местами положения задач.
  • Сохраняйте естественный линейный порядок повествования в ноутбуке сверху-вниз.
  • Отвечайте на вопросы, а также добавляйте новые ячейки в предложенных местах, которые обозначены <...>.
  • В markdown-ячейка, содержащих описание задачи, находятся специальные отметки, которые запрещается модифицировать.
  • При нарушении данных правил работа может получить 0 баллов.

Перед выполнением задания посмотрите презентацию по выполнению и оформлению домашних заданий с занятия 2.

Баллы за задание:

  • Задача 1 — 40 баллов
  • Задача 2 — 70 баллов
  • Задача 3 — 20 баллов

In [ ]:
# Bot check

# HW_ID: phds_hw3
# Бот проверит этот ID и предупредит, если случайно сдать что-то не то

# Status: not final
# Перед отправкой в финальном решении удали "not" в строчке выше
# Так бот проверит, что ты отправляешь финальную версию, а не промежуточную
In [ ]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression, LinearRegression
from sklearn.metrics import accuracy_score, balanced_accuracy_score

from typing import List, Callable, Tuple

import seaborn as sns
sns.set_theme('notebook', font_scale=1.2, palette='Set2')

from tqdm.notebook import tqdm

Задача 1¶

Рассмотрим модель логистической регрессии. Признаки объекта представимы в виде d-мерного вектора x∈Rd, класс имеет бернуллиевское распределение Y∼Bern(μθ(x)). Мы делаем следующее предположение о зависимости параметра вероятности от признаков μθ(x)=σ(xTθ)=11+e−xTθ.

При добавлении регуляризации к модели логистической регрессии оптимизируемый функционал принимает вид

F(θ)=−n∑i=1[Yilogσ(θTxi)+(1−Yi)log(1−σ(θTxi))]+λθTθ

  1. Выпишите формулы градиентного спуска (GD) и стохастического градиентного спуска (SGD).

  2. Покажите, что F(θ) — выпуклая функция по θ и, как следствие, имеет единственный экстремум, являющийся глобальным максимумом. Указание. Посчитайте гессиан (матрицу вторых производных) и покажите, что она положительно определена.

  3. Опишите, как может вести себя решение при отсутствии регуляризации, то есть при λ=0

Задача 2¶

Введение¶

Вы когда-нибудь задумывались, что такое запах и вкус? Как именно мы различаем запахи? Ведь это даже не физическая величина: наш нос не измеряет концентрацию примесей в воздухе, их парциальное давление или что-то ещё. Запах — есть субъективная реакция мозга на химическое взаимодействие молекул с рецепторами в слизистой оболочке носа.

Рассмотрим, на первый взгляд, забавное устройство — электронный нос (e-nose). Этот "прибор", буквально, повторяет принцип работы человеческого носа и состоит из трёх основных компонентов:

  1. Массив датчиков, чувствительных к наличию различных примесей в газе
  2. Канал доставки, по которому газ поступает к датчикам
  3. Система обработки данных, анализирующая сигналы от датчиков и сопоставляющая их с известными шаблонами. Если запах уже "знаком" системе, она может его распознать

image.png

Картинка из Википедии

Спектр возможных применений подобных устройств может быть весьма широк: от контроля качества скоропортящихся продуктов до обнаружения взрывчатых веществ, вирусов, бактерий или токсичных газов.

Сами сенсоры бывают разными. Наиболее современные разработки основаны на технологиях с применением графена и углеродных нанотрубок. Более доступной, но достаточно эффективной альтернативой являются полупроводниковые сенсоры на основе оксидов металлов. Их принцип работы основан на изменении проводимости при адсорбции молекул газа на поверхности датчика. Измеряя сопротивление таких сенсоров, можно определить состав и даже концентрацию примесей в воздухе.

Проблема¶

Одной из самых сложных проблем в этой области считается дрейф сигнала датчиков. Он возникает из-за физико-химического воздействия окружающей среды: со временем на сенсоры оседают посторонние молекулы, изменяя структуру и свойства материала. Системе обработки становится всё сложнее анализировать данные, поскольку она откалибрована (обучена) на свежих сенсорах, без учёта дрейфа. Это, очевидно, приводит к постепенному ухудшению точности измерений.

Для изучения этого явления был проведён долгосрочный эксперимент, в ходе которого на протяжении трёх лет фиксировались "отпечатки запахов" шести летучих органических соединений при различных концентрациях: аммиака, ацетальдегида, ацетона, этилена, этанола и толуола. Подробное описание методики, сбора данных и выделения признаков можно найти в оригинальной статье, а также на странице с данными. Авторы исследования разработали и показали эффективность нового (на момент выхода статьи) способа учёта дрейфа в классификаторах газов.

Данные¶

Исследуемый "нос" состоял из 16 сенсоров. Для каждого сенсора снималась зависимость изменения сопротивления от времени.

  1. Максимальное изменение сопротивления в ходе одного измерения (абсолютное и относительное) — 2 признака
  2. Максимум/минимум экспоненциального скользящего среднего (EMA) при разных параметрах сглаживания — 6 признаков
Подробнее ✍️

Для последовательности временных отсчётов r1, r2, r3... экспоненциальное скользящее среднее рассчитывается как yk=(1−α)yk−1+α(rk−rk−1), Параметр α∈(0,1) определяет степень "забывания" старых значений: чем он меньше, тем сильнее сглаживается сигнал. В нашем случае EMA рассчитывалось отдельно для возрастающего и убывающего участков сигнала при α=10−1,10−2,10−3. В качестве признаков взяты, соответственно, максимальное и минимальное значения за время измерения. Подробнее см. в статье.


Итого, имеется 8⋅16=128 признаков

1. Загрузка и подготовка данных¶

1.1 Смотрим на данные¶

Скачайте датасет. Если открыть любой файл в текстовом редакторе, вы увидите, что данные хранятся в разреженном формате svmlight. Для его загрузки в sklearn предусмотрена специальная функция. Воспользуйтесь ей:

In [ ]:
from sklearn.datasets import load_svmlight_files

file_paths = [f'batch{n}.dat' for n in range(1, 11)]

data = <...>

type(data), len(data), type(data[0]), type(data[1])
Out[ ]:
(list, 20, scipy.sparse._csr.csr_matrix, numpy.ndarray)

Вы получили список из 20 элементов, в котором чередуются матрицы признаков и векторы таргета. Матрицы признаков имеют тип разреженной матрицы. Преобразуем их в pandas.DataFrame и присвоим осмысленные имена столбцам:

In [ ]:
# Формируем список имён признаков
feature_names = [
    f"{sensor}_{suffix}"
    for sensor in range(1, 17)
    for suffix in ["dR", "norm_dR"] +
                  [f"max_ema_{alpha}" for alpha in [0.1, 0.01, 0.001]] +
                  [f"min_ema_{alpha}" for alpha in [0.1, 0.01, 0.001]]
]

# Разделяем данные: X (признаки) и y (таргет)
X_raw, y_list = data[::2], data[1::2]

# Преобразуем разреженные матрицы в DataFrame
X_list = [pd.DataFrame(X.toarray(), columns=feature_names) for X in X_raw]

X_list[0].head()

Замечание. Первое число в названии означает номер сенсора.

Почему данные представлены в таком виде? Согласно описанию эксперимента, количество измерений для разных газов варьировалось от месяца к месяцу. В некоторые периоды отдельные газы вовсе не использовались. Чтобы сбалансировать распределение классов, данные сгруппировали в 10 батчей, где разные газы представлены более равномерно.

Проверим, насколько хорошо это удалось: напишите функцию, которая визуализирует распределение классов в каждом батче, и примените её к данным.

Подсказка ✍️

Удобно воспользоваться одной из функций: sns.countplot, sns.barplot, plt.bar


In [ ]:
def visualize_classes(target_batches: List[np.ndarray]):
    """Визуализирует распределение классов во всех батчах

    Args:
        target_batches (List[np.ndarray]): список таргетов
    """
    
    <...>


visualize_classes(y_list)

Вывод:

Выведите количество экспериментов в каждом из батчей.

In [ ]:
 

Вывод:

1.2 Учитываем время¶

Поскольку объекты нашей выборки упорядочены по времени, имеет смысл добавить признак experiment_id, описывающий номер эксперимента. Причём, нумерация должна быть сквозной: если последний эксперимент 1-го батча имеет номер N, то первый эксперимент второго батча должен иметь номер N+1. Создайте этот столбец и поместите его в начало каждой таблицы с признаками:

In [ ]:
 

Проверьте, всё ли правильно вы сделали:

In [ ]:
# Создаем списки для первой и последней строки таблицы
first_ids = [batch['experiment_id'].iloc[0] for batch in X_list]
last_ids = [batch['experiment_id'].iloc[-1] for batch in X_list]

# Создаем DataFrame для отображения
result = pd.DataFrame(
    [first_ids, last_ids],
    index=['First experiment_id', 'Last experiment_id'],
    columns=[f'batch {i+1}' for i in range(len(X_list))]
)

result

1.3 Определяемя с таргетом¶

Задачу многоклассовой классификации решать сложно. Но можно свести её к бинарному случаю. Для этого выберите один газ и обозначьте его как класс 1. Все остальные газы отнесите к классу 0.

📌 Замечание. Вместо 0 и 1 можете придумать более понятные названия.

In [ ]:
gas_id = <...>

<...>

Такая стратегия имеет название "One-VS-All". При помощи написанной вами ранее функции выведите распределение новых классов по батчам:

In [ ]:
 

Далее, при построении классификатора, учитывайте диcбаланс: используйте взвешенную версию логистической регрессии (см. параметр class_weight) и метрику balanced_accuracy_score для оценки качества модели. ❗Это очень важно❗

1.4 Что собираемся делать?¶

Основное назначение этого датасета — разработка алгорима машинного обучения, способного учиывать дрейф датчиков. Авторы исследования предложили решение, основанное на построении взвешенного ансамбля линейных классификаторов.

Мы же пока сосредоточимся на более простых вещах. В качестве обучающей выборки будем использовать первый батч, и посмотрим, как будет деградировать точность логистической регрессии с течением времени. Выделите обучающий батч и список тестовых батчей.

In [ ]:
 

1.5 Стандартизация¶

Обучите StandardScaler на обучающем батче и преобразуйте его. Преобразуйте тестовые батчи.

In [ ]:
 

2. Есть ли мультиколлинеарность?¶

Хорошо ли обусловлена обучающая матрица признаков? О чём это говорит?

In [ ]:
 

Вывод:

Для наглядности, выведите матрицу корреляций. Что означают элементы матрицы? А в нашем случае?

In [ ]:
 

Ответ:

3. Модельки, модельки, модельки...¶

3.1 Самый популярный класс¶

Найдите самый популярный класс в обучающем батче и посчитайте точность ответа на тестовых данных только этим классом — константой. Отличается ли взвешенная точность от обычной?

In [ ]:
 

Вывод:

3.2 Модель логистической регрессии без регуляризации¶

Реализуем функции для обучения, тестирования и отображения деградации метрики качества ответов модели со временем (заполните пропуски):

In [ ]:
def visualize_metric_fall(scores: List[float]):
    """Визуализирует метрику на обучающем и тестовых батчах с пострением
    прямой характерного спада

    Args:
        scores (List[float]): значения метрики для различных батчей
    """

    n = np.arange(len(scores) - 1)

    linreg = LinearRegression(fit_intercept=True)
    linreg.fit(n.reshape(-1, 1), scores[1:])

    a, b = linreg.coef_[0], linreg.intercept_

    plt.figure(figsize=(8, 5), tight_layout=True)

    plt.plot(n, n * a + b, label=f'${a:.2f}n + {b:.2f}$')
    plt.scatter(n, scores[1:], label="Тестовые батчи")
    plt.scatter(-1, scores[0], label='Обучающий батч')
    plt.xlabel('Номер тестового батча')
    plt.ylabel('Значение метрики')
    plt.legend()


def evaluate_accuracy_decay(
    model: LogisticRegression,
    X_train: pd.DataFrame,
    y_train: np.ndarray,
    X_test_batches: List[pd.DataFrame],
    y_test_batches: List[np.ndarray],
    metric: Callable[[np.ndarray, np.ndarray], float],
    visualize: bool = False
) -> List[float]:
    """Обучает модель на X_train и y_train, тестирует на X_test_batches.
    Возвращает список значений метрики на каждом тестовом батче.

    Args:
        model (BaseEstimator): модель логистической регрессии sklearn
        X_train (pd.DataFrame): обучающий батч (признаки)
        y_train (np.ndarray): обучающий батч (таргет)
        X_test_batches (List[pd.DataFrame]): тестовые батчи (признаки)
        y_test_batches (List[np.ndarray]): тестовые батчи (таргет)
        metric (Callable[[np.ndarray, np.ndarray], float]): метрика для оценки качества
        visualize (bool): если True, визуализировать изменения метрики

    Returns:
        List[float]: список значений метрики на тестовых батчах
    """
    
    # Обучение модели
    <...>

    # На первом месте будет стоять значение метрики на обучающей выборке
    scores = [metric(y_train, model.predict(X_train))]

    # Вычисляем метрику на каждом тестовом батче и добавляем в scores
    <...>
    

    if visualize:
        visualize_metric_fall(scores)

    return scores

Обучите классическую логистическую регрессию без регуляризации и визуализируйте точность ответа на трейне и тесте. Свободный коэффициент необходимо исключить из модели.

Подсказка ✍️

Чему равен аргумент penalty по умолчанию?


In [ ]:
 

Вывод:

Далее выполняйте тестирование только на первом тестовом батче.

3.2 Логистическая регрессия с регуляризацией¶

За что отвечает гиперпараметр C у класса LogisticRegression?

Ответ:

Вам необходимо исследовать зависимость от C следующих величин:

  1. Accuracy на трейне
  2. Accuracy на тесте
  3. Коэффициенты модели

Чтобы не приходилось постоянно обучать модели при одних и тех же сетках C, предлагается написать функцию, которая будет принимать на вход вид штрафа penalty, границы диапазона C, и саму выборку. На каждой итерации вычисляйте все величины и сохраняйте в виде списков. Для мониторинга времени работы используйте функцию tqdm. Пример использования:

from tqdm.notebook import tqdm
for C in tqdm(C_grid):
    <тело цикла>

Не забудьте также про имеющийся дисбаланс классов.

In [ ]:
def train_alpha_grid(
    min_log_C: float,
    max_log_C: float,
    resolution: int,
    X_train: pd.DataFrame,
    y_train: np.ndarray,
    X_test: pd.DataFrame,
    y_test: np.ndarray,
    penalty: str,
    solver: str = 'newton-cholesky',
    max_iter: int = 100
) -> Tuple[np.ndarray, List[List[float]], List[float], List[float]]:
    """Обучает модель LogisticRegression для разных значений параметра регуляризации C,
    сохраняет коэффициенты, вычисляет accuracy на обучающей и тестовой выборках.

    Args:
        min_log_C (float): минимальное значение log10(C) для сетки.
        max_log_C (float): максимальное значение log10(C) для сетки.
        resolution (int): число точек на сетке C.
        X_train (pd.DataFrame): обучающая выборка (признаки).
        y_train (np.ndarray): отклик на обучающей выборке.
        X_test (pd.DataFrame): тестовая выборка (признаки).
        y_test (np.ndarray): отклик на тестовой выборке.
        penalty (str): тип регуляризации ('l1', 'l2', 'elasticnet', 'none').
        solver (str, optional): метод оптимизации параметров модели. По-умолчанию 'newton-cholesky'.
        max_iter (int, optional): максимальное количество итераций для оптимизации. По-умолчанию 100.

    Returns:
        Tuple[np.ndarray, List[List[float]], List[float], List[float]]:
            - C_grid (np.ndarray): сетка значений C,
            - coefs_list (List[List[float]]): список коэффициентов для каждого значения C,
            - baccuracy_train_list (List[float]): список balanced accuracy на обучающей выборке для каждого значения C,
            - baccuracy_test_list (List[float]): список balanced accuracy на тестовой выборке для каждого значения C.
    """

    <...>

Проведите эксперимент для 3-х разных моделей логистической регрессии с различными типами регуляризации:

  1. L1-регуляризация
  2. L2-регуляризация
  3. Комбинированная регуляризация с параметром l1_ratio=0.5.

Рекомендации

  • Подберите диапазоны значений для гиперпараметра C. Не берите слишком узкие, чтобы видеть на графике всю картину. Для слишком широких границ придётся брать больше точек.
  • Вам не нужна очень частая сетка гиперпараметра C. При отладке кода можно вообще использовать сетку из 2-3 значений.
  • Вы можете столкнуться с различными ошибками и warning-ами (например, неверный solver, отсутствие сходимости, и т.д.) Постарайтесь настроить гиперпараметры модели таким образом, чтобы ошибки исчезли, а количество предупреждений было минимальным.
  • Для ускорения работы программы можно использовать параллельные вычисления, передав аргумент n_jobs в модель. Это особенно полезно при запуске на локальном компьютере, так как Google Colab предоставляет лишь два ядра ЦПУ по-умолчанию.
In [ ]:
 

Нарисуйте треки коэффициентов моделей в зависимости от C. Легенду можно сделать общую, если все графики помещаются на экране. Отразите в ней наименования признаков для соответствующих коэффициентов. Сделать красиво могут помочь заметки отсюда.

In [ ]:
 

Вывод:

Нарисуйте зависимости точности предсказания от C на обучающей и тестовой выборках. Скомпонуйте всё на 2-3 графиках. Горизонтальными линиями, отметьте точность модели без регуляризации на трейне и тесте.

In [ ]:
 

Вывод:

Установите гиперпараметры логистической регрессии, при которой достиглось оптимальное качество. Визуализируйте точность её работы в зависимости от номера батча. Сравните с моделью без регуляризации. Сделайте выводы.

📌 Замечание. Имейте в виду, что точность на обучающей выборке должна оставаться близкой к 1. Это несколько противоречит классическому подходу, когда для борьбы с переобучением стараются одновременно улучшать метрику на тесте и ухудшать на трейне. В нашем случае следует предполагать (и надеяться), что система электронного носа спроектирована правильно, и газы хорошо отделяются друг от друга, пока датчики ещё "свежие". В таком случае три кита хорошей модели это:

  1. Точность на трейне близка к 1
  2. Высокая точность на первых тестовых батчах
  3. Слабый спад точности с течением времени.
In [ ]:
 

Сделайте общий вывод по задаче. Есть ли пользва от регуляризации с точки зрения метрики? Получается ли удовлетворить всем трём перечисленным условиям одновременно?

📌 Замечание. Учитывайте, что мы свели задачу классификации 6 газов к бинарному случаю. Поэтому результаты могут сильно различаться при разных разбиениях. Это вполне логично: дрейф сенсоров крайне плохо предсказуем, и не обязан одинаково влиять на способность системы различать разные газы.

Вывод:

Задача 3 (продолжение)¶

Продолжаем работать с этим датасетом.

1. Число обусловленности¶

Исследуйте зависимость числа обусловленности от параметра C для L1-регуляризации. Постройте соответствующий график.

In [ ]:
 

Вывод:

2. Предсказание вероятностей¶

Исследуйте распределение предсказываемых вероятностей для логистической регерессии с регуляризацией и без.

Для начала, реализуйте функцию, которая будет обучать логистическую регрессию с наиболее оптимальным C и возвращать предсказание вероятности для касса 0. Для этого у всех классификаторов sklearn предусмотрен метод predict_proba.

In [ ]:
def get_proba_distr(
    X_train: pd.DataFrame,
    y_train: np.ndarray,
    X_test: pd.DataFrame,
    C: float = np.inf
) -> Tuple[np.ndarray, float]:

    <...>

Сравните гистограммы распределения предсказанных вероятностей для двух указанных моделей. Используйте первый тестовый батч.

In [ ]:
 

Вывод:

Постройте аналогичные гистограммы для ещё двух тетовых батчей. Сделайте выводы. Меняет ли регуляризация распределение вероятностей?

In [ ]:
 

Вывод: